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Abstract Most soft tissues possess an oriented architecture of collagen fiber bun- 
dles, conferring both anisotropy and nonlinearity to their elastic behavior. Trans- 
verse isotropy has often been assumed for a subset of these tissues that have a sin- 
gle macroscopically-identifiable preferred fiber direction. Micro- structural studies, 
however, suggest that, in some tissues, collagen fibers are approximately normally 
distributed about a mean preferred fiber direction. Structural constitutive equations 
that account for this dispersion of fibers have been shown to capture the mechani- 
cal complexity of these tissues quite well. Such descriptions, however, are compu- 
tationally cumbersome for two-dimensional (2D) fiber distributions, let alone for 
fully three-dimensional (3D) fiber populations. In this paper, we develop a new 
constitutive law for such tissues, based on a novel invariant theory for dispersed 
transverse isotropy. The invariant theory is based on a novel closed-form ‘splay in- 
variant’ that can easily handle 3D fiber populations, and that only requires a single 
parameter in the 2D case. The model is polyconvex and fits biaxial data for aor- 
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tic valve tissue as accurately as the standard structural model. Modification of the 
fiber stress-strain law requires no re-formulation of the constitutive tangent matrix, 
making the model flexible for different types of soft tissues. Most importantly, the 
model is computationally expedient in a finite-element analysis. 

Key words nonlinear elasticity - invariant theory - anisotropy 


1 Introduction 

Constitutive modeling of soft tissues is an important prerequisite for the com- 
putational analysis of the systems-level response of biological structures to ap- 
plied loads and pressures. Computational modeling in many cases is the only way, 
short of implantation, to evaluate the mechanical response of tissues to the full 
3D mechanical environment seen in vivo. Moreover, computational modeling per- 
mits the analyst to directly query the relationship between structure and stress, and 
to explore the design space of newly proposed, implants, bioprostheses or tissue- 
engineered constructs. Numerical models of tissue systems can be computationally 
intensive. This is particularly true of tissues that require complex material models, 
and systems that involve fluid-structure interactions. The material models used in 
these systems must not only be accurate, they must also be computationally effi- 
cient. 

The past few years have seen an increased interest in nonlinear continuum me- 
chanics as a framework for describing the mechanical behavior of soft tissues. The 
now well-established mathematics of this field provides a perspective from which 
rigorous, thermodynamically-reasonable constitutive equations can be proposed — 
a characteristic that was lacking in many earlier ad hoc tissue descriptions. The 
geometric and material nonlinearities seen in tissues fit well into this framework, 
and anisotropy is readily handled by the theory of invariants [Spencer, 1972]. In 
addition, constitutive equations posed within this framework can call upon de- 
veloped computational techniques, permitting the exploration of tissue-level and 
organ-level mechanics involving finite deformations. An excellent reference for 
both nonlinear continuum mechanics and related variational principles, as they 
apply to soft tissues, is the textbook by [Holzapfel, 2000] 

Structural constitutive equations that view soft tissues as statistically-oriented 
distributions of fibers come from a different tradition based on the experimen- 
tal observation of fiber dispersion, or splay. For example, [Sacks, 2003] showed 
a correlation between the intensity distribution of the small-angle light scatter- 
ing of bovine pericardium and both fiber orientation and dispersion. Similarly, 
[Holzapfel et al., 2002] showed that smooth muscle cells in arterial media are sta- 
tistically oriented around two opposing helical directions. From these cellular ori- 
entations, collagen fiber orientations were inferred. An earlier observation of fiber 
dispersal in bovine pericardium was documented by [Zioupos and Barbenel, 1994]. 
The need to account for fiber dispersion architectures when modeling soft tissues 
was first addressed by [Lanir, 1983]. Since then, structural models with disper- 
sion have been proposed for passive myocardium [Horowitz et al., 1988], lung 
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[Mijailovich et al., 1993], heart valves [Sacks, 2000], aorta [Wuytsetal., 1995], 
and tendon and ligament [Hurschler et al., 1997]. To model fiber splay, typically 
either a von Mises [Hurschler et al, 1997] or Gaussian [Sacks, 2000] distribution 
is adopted. 

Mathematically, the inclusion of such probability distribution functions into 
constitutive models presents little challenge, but computationally, they are cum- 
bersome. The probability density function acts as a weighting function for a fiber 
stress-strain rule. Functionally, in the 2D case, this can be represented as 

*h 

S = J S f (6) R(6) N(0) ® N(0) dd, (1) 

~ r h 

where S is the second Piola-Kirchhoff stress, S f (6) is a fiber stress-strain rule, 
R(6) is a probability distribution function, and N(0) ® N(0) governs the orientation 
of a fiber family in the 2D plane. What is important to note is that the product 
of fiber stress and the weighting function must be integrated over a full semi- 
circle in the 2D case, Eq. (1), or an equivalent hemisphere in the 3D case. These 
integrals have traditionally been evaluated numerically. Numerical experiments on 
one such equation for mitral- valve tissue suggested that eighteen discrete intervals 
were necessary to capture the full range of constitutive behavior [Einstein, 2002] . 
This is mechanically equivalent to having eighteen weighted fibers in the plane. 
A single integral of this kind in the constitutive equation implies two integrals 
to be evaluated in the constitutive tangent matrix. Likewise, 3D fiber dispersion 
would require two integrals to be evaluated for stress, and four for the constitutive 
tangent matrix. In the computational model, these operations are evaluated at every 
iteration, of every time step, at every Gauss point, for every finite element, in some 
geometric model of interest. 

In this paper, we develop an alternative constitutive law for tissues whose col- 
lagen fiber populations are statistically distributed. To address the problem of com- 
putational cost, the integral in Eq. (1) is replaced with a novel closed-form ‘splay 
invariant’ that requires a single parameter in the 2D case, and a single operation 
per iteration. To evaluate the model, we compare its correlative capability against 
biaxial data for aortic-valve tissue. In addition, we compare the capabilities of 
our model against those of a published structural model [Billiar and Sacks, 2000b, 
Einstein, 2002] that is based on the paradigm in Eq. (1) and which fits the data 
quite well. 


2 Theory 

Consider a mass point originally given by the set of coordinates X = (Xf, X 2 , X3) 
assigned at an arbitrary reference time of t 0 . At current time t, this mass element is 
located by a different set of coordinates x = (jti , jc 2 , * 3 ). Let the motion of this mass 
point through space be described by a one-parameter family (in time) of locations 
considered to be continuous and sufficiently differentiable to allow the definition 
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of a deformation gradient 
3xi 

Fij(to, t ) = — or in tensor notation, F(fo, t) = \dx r ldX c \, (2) 
aXj 

where indices r and c denote row and column, and have values 1, 2, 3. The ability 
to invert this field (i.e., F _1 (fo,0 = \dX r ldx c \ such that F -1 F = FF -1 = I, where 
I = is the identity tensor and 6 t j denotes the Kronecker delta) ensures that a 
given particle cannot occupy two locations at the same instant in time, and that two 
discrete particles cannot occupy a single location at any given moment in time. 

Affiliated with the deformation gradient defined in Eq. (2) are the left- and 
right-deformation tensors 1 defined by 

B = FF t and C = F T F, (3) 

respectively, where ‘ T ’ implies matrix transpose (e.g., B(x; to, t) = |F r jfcF c fc| and 
C(X; to, t ) = iFkrF^ with the repeated index k being summed from 1 to 3 in 
the usual manner). Inverses B -1 and C" 1 exist because the tensor fields B and C 
are symmetric positive-definite. The left-deformation tensor B appears in Eulerian 
constructions, while the right-deformation tensor C appears in Lagrangian con- 
structions. 


2.1 Invariants 

Consider a vector ao = |a r (X; ?o)J of length ao- ao = a*(X; to) a*(X; to) = 1 that lies 
tangent to a material line of strength (e.g., a fiber) in the reference state to- After 
a deformation, this material line will have stretched by an amount T(/o, 0 that is 
quantified by A 2 = ao - Cao = a*(X; to) Qr(X; to, 0 a^(X; to). As such, there exists a 
spatial vector a = [a r (x; ?)]] with L2 norm ||a ||2 = (a ■ a) 1/2 = A that in the deformed 
state t of Green’s metric C relates to material vector ao via the mapping a = Fao, 
which is the transformation law of a polar vector. 2 

From the classic theory of invariants [Spencer, 1972], there are five invariants 
that are needed to describe a material with transverse isotropy (i.e., a material with 
a single fiber family); they are: 

// = tr C, I n = i((tr C) 2 - tr(C 2 )), / m = det C = (det F) 2 , (4) 

I iv - ao • Cao = a • a = A 2 , Iv = ao • C 2 ao = a ■ Ba, (5) 

1 Tensors B = FF T = V 2 and C = F T F = U 2 are often referred to as the left and right 
Cauchy-Green deformation tensors, respectively, because V and U are called the left- and 
right-stretch tensors, so named because of the polar decomposition F = VR = RU wherein 
RR t = l. We prefer to call B and C the left- and right-deformation tensors. Historically, 
[Cauchy, 1827, pp. 60-69] used B _1 (sometimes expressed as c) and [Green, 1841] used C, 
while [Finger, 1894] introduced their duals, B (sometimes expressed as b) and C" 1 . There- 
fore, naming these fields after Cauchy and Green seems an injustice to Finger, especially 
since Finger introduced B into the literature. 

2 The transformation law b = F“ T b 0 applies to axial vectors whose norms ||b || 2 = (b-b ) 1/2 
are measured against Finger’s metric C -1 . We have no need for axial vectors in this work. 
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where trC = Ckk is the trace of C, and detC = j{trC 3 - trC[(trC) 2 - 3 tr C 2 ]} 
denotes its determinant. The three invariants in Eq. (4) account for isotropic effects, 
while the two invariants in Eq. (5) account for anisotropic effects. 

Invariants Ijy and ly idealize all fibers in a family as being parallel, which is 
not indicative of soft tissues. The fiber architectures of soft tissues are splayed. 
We therefore seek an alternative pair of invariants — call them and 7<y> — that 

can replace Ijv and ly for dispersed fiber architectures, yet analytically reduce to 
Ijy and ly in the absence of fiber dispersion. It is sufficient to define these new 
invariants as 


I (iv) = tr(FKF T ) - tr(KC) and /<y> = tr(CKC), (6) 

where K(X; r 0 ) is constrained such that K a 0 ® a 0 = ||a r (X; to)a c (X\ f 0 )J in the 
absence of splay. Tensor K is a material constant. The main objective of this paper 
is to derive such a K appropriate for describing the anisotropy caused by fiber 
dispersion. 


2.2 Elasticity 


The strain-energy density per unit mass, when written in the Lagrangian frame, is 
given by [Lodge, 1974, pp. 194-195] 


dW = -^-tr(SdC), 

2po 


(7) 


where dW(X\to,t, dt) represents the work done on a material element of mass den- 
sity g = p(x;f), with go = £>(X;to)- Work is caused by an imposed displacement 
acting on the mass element, manifested here as the strain increment ^dC(X; to, t, dt). 
The material responds to this displacement through the creation of forces, thereby 
producing a state of stress S(X; to, t). 

The second Piola-Kirchhoff stress tensor S pulls forward into the Eulerian 
frame [Holzapfel, 2000, pp. 82-84] becoming the Cauchy stress tensor T(x; t) via 
the well-known mapping S = §F -1 TF -t . Formula ^ = detF follows from the 
conservation of mass. Green strain E(X; to, 0, defined by E = i(C-l), has an incre- 
mental change dE(X; to, t, dt) of dE = -dC = F T EF, wherein E(x; t, t+dt) = ^(0—1) 
with C(x; t, t+dt) = F T F given that F(t, t+dt) = \dx r {t+dt)ldx c {t)\. 

An elastic solid is defined by the constitutive law [Leonov, 2000] 


S = 2qq 


dW(T, C) 
dC 


( 8 ) 


with an isochoric constraint of det F = 1 also applying whenever the material is 
incompressible. Thermodynamics requires W to be a function of both temperature 
T and deformation C. Because the human body maintains a nearly isothermal state, 
the temperature dependence of tissues is usually neglected in their analysis. The 
strain-energy density W will also depend on any number of material constants that 
may appear as scalar, vector, or tensor fields. 
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Adopting the invariants in Eqs. (4 & 6) as our integrity basis, the general consti- 
tutive equation for a transversely isotropic elastic solid with splay, when expressed 
in the Lagrangian frame, yields the constitutive equation 3 

S = ?Wj + Ii<Wji)\ - + ImWin C~ l + W { w } K + W <V >(CK + KC), (9) 

that, when pushed forward into the Eulerian frame, becomes 

r = (<W tI + IfWj) B - W n B 2 + l m <W M I + W (IV ) A + ^V <V) (BA + AB), (10) 

where *W/ = 2godW/dIj, etc., and A(x;to, 0 = FKF T . Tensor r(x; to, t) = is 
known as the Kirchhoff stress. 

It is a straightforward matter to extend any, existing, transversely isotropic, 
integrity basis into an equivalent basis that accounts for splay by applying an ap- 
propriate mapping. For example, the classic integrity basis {//, In, Im, In, Iv) maps 
into the basis \Pi,Pu,Pm,Piv,Pv } derived by [Criscione et al., 2001] via their for- 
mulae (5.3 & 5.4). To introduce splay into their integrity basis, one simply applies 
these same mappings, but with {//, In, Im, I(iv), /<v>} replacing {//, l n , Im, In, Iv}, 
which produces {J3i,fi(n),f3(m),fi{N),/3(V)} as its equivalent set of splay invariants. 

The challenges that lie ahead are: i ) to establish a tensor field K which is ap- 
propriate for splayed fiber architectures, and ii) to arrive at a simple set of five 
scalar-valued gradients {' W I , < W y n, r W t m,W(N), r W < (v >} that allows Eqs. (9 & 10) 
to aptly describe some known set of experimental data. 


3 Anisotropic Stiffness 

For purposes of assessing the effect of fiber orientation on stiffness, it is useful to 
switch from the global coordinate system (X, Y, Z) with base vectors {ex, ey, e z } to 
a local or intrinsic coordinate system (1, 2, 3) with base vectors {ei, e2, 63}. These 
local coordinates are selected so that the unit vector in the 1 -direction (i.e., ei) is 
coaxial with the mean direction of fiber orientation a 0 , while the unit vectors in 
the 2- and 3-directions lie in the transverse plane. Because both the intrinsic and 
global coordinate systems are considered to be rectangular Cartesian, there exists 
an unique, orthogonal, rotation matrix Q such that 


ex = Qei, = Qe 2 , e z = Qe 3 , 


(11) 


where Q T Q = I with detQ = 1. Matrix Q represents a rigid-body rotation. 


3 The following well-known tensor derivatives are useful in the derivation of constitutive 
formulae: 


dtrZ 3trZ -1 

~dZ~ ~ dZ 


-z- T z- T , 


ddetZ 

dZ 


= det(Z)Z _T . 
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3.1 Two-Dimensional Splay Invariants 


For transversely isotropic 2D membranes, the unit vector e f in the intrinsic coor- 
dinate frame that locates a specific fiber tangent within a given fiber distribution is 
described via polar coordinates 


e f (0) = cos(0)ei + sin(60e 2 , (12) 

where angle 6 orients e/ with respect to the mean-fiber direction ei (i.e., to ao) in 
the plane of a membrane (viz., the 12 plane). 

We adopt the Gaussian formulation primarily because it allows analytic solu- 
tions (cf. Apps. A & B). The cosine of the angle between the mean-fiber direction 
ei and that of an individual fiber e/ is given by the inner product eye/, which from 
Eq. (12) is just cos#. The Gaussian distribution governing a single fiber family 
(i.e., transverse isotropy), when expressed in the local coordinate system (1,2, 3) 
with base vectors {ei, e 2 , e 3 }, is therefore simply 

^ exp (5)’ 

with a being a standard deviation in the angle of fiber dispersion about the mean- 
fiber direction ei. Because the local coordinate direction ei is coincident with the 
mean-fiber direction ao, there is no angle between ei and ao, and as such, the mean 
angle for this distribution, which is usually denoted by p, is identically zero — a 
direct consequence of selecting the intrinsic coordinate system that we did. 

In the spirit of [Lanir, 1983], we propose the following definition for our forth 
invariant for 2D splay: 


v> 


17 1 2 



exp(-fl 2 / 2 < 7 2 ) 

yl2ntri{nl2vj2g) 


e/(0) -Q t CQ e/(0)d0, 


(13) 


where g is akin to a, except that g is not a standard deviation whose units are radi- 
ans. A like expression with C 2 replacing C defines /<y>. The error function erf(jc) is 
introduced into these formulas for reasons that will be made clear shortly. These in- 
variants satisfy the required limits: lim^ 0 I{iv)(g) = Ijv and lim^o l {V )(g) = h- In 
Eq. (13), the right-deformation tensor C is rotated into the local coordinate frame 
(1, 2, 3) by the mapping Q T CQ, or equivalently, the fiber tangent vector e/ is ro- 
tated into the global coordinate frame ( X , Y, Z) by the mapping Q e/, in accordance 
with Eq. (11). 


3.2 Three-Dimensional Splay Invariants 

For transversely isotropic 3D tissues, e/ is located via spherical coordinates 

<p) = cos(0) ei + sin(0)(cos(0) e 2 + sin(0) e 3 ), (14) 

based on the geometry presented in Fig. 1 . Here angles 0 and 0 orient e / with 
respect to the mean-fiber direction ei of the embedded frame. 
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Figure 1 The diagram on the left relates the global coordinates ( X , Y, Z ) to the local coor- 
dinates (1,2, 3), selected so that the mean-fiber direction ^ in the Lagrangian frame aligns 
with the 1 axis. The diagrams on the right illustrate how the unit vector e/ for a specific fiber 
within a fiber distribution of a 3D tissue is oriented with respect to the mean-fiber direction 
a 0 via angles 6 and 


For 3D splay with transverse isotropy, our forth invariant is defined by 
exp(- 0 2 / 2 < 7 2 ) 


'W 


= JJ 

0 -*/2 


qiffnff 2 txi{n I 


e f {6, 0) • Q CQ e j{9, <p) d# d0. (15) 


A like formula with C 2 replacing C defines /<y>. The Gaussian distribution is in- 
dependent of (f) in these invariants, because of an assumption of isotropy in the 
transverse plane. Once again, lim^o /(/vote) = Ijv and lim f _> 0 /<v>te) = V 

The fact that a right circular cone is used to describe the angular dispersion of 
individual fibers implies a radial symmetry in the transverse plane, which is con- 
sistent with the notion of transverse isotropy. An elliptic cone could be employed 
if radial symmetry were deemed inappropriate. We will discuss this case later, but 
we do not derive it. 


3.3 Intrinsic Anisotropic Stiffness 


We postulate the existence of a material-constant tensor field that we denote as 
k, which serves as a relative (i.e., normalized) stiffness matrix associated with the 
anisotropic facets of material geometry. For 2D splay, this field is given by 


"h 

*( f ) = J 


~ x l2 


exp(-# 2 /2^ 2 ) 
g yfln exi{n /lyjlg) 


effO) ® e/(0)d 0, 


(16) 
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and for 3D splay with transverse isotropy, it is given by 


2 n ” h 


*(?) = J J 


0 ~«h 


exp(-0 2 /2^ 2 ) 
s(2n) 3/l erf(^’/2'\/5^) 


ey (9, 0) <s> ey (6, 0) d 6 d 0, 


(17) 


where the global stiffness K introduced in Eq. (6) relates to this intrinsic stiffness 
k via 

K = QkQ t . (18) 

Both of these k matrices obey the required property: lim s ~_,o Q k(<t)Q t = ao ® ao, 
and both of the k matrices are symmetric; therefore, their affiliated K matrices are 
symmetric, too. Analytic solutions to Eqs. (16 & 17) are provided in Apps. A & B, 
respectively. The ability to analytically solve for the anisotropic stiffness K means 
that this theory will be efficient when implemented into finite-element codes. In 
contrast, the models of[Hurschler et al., 1997], [Sacks, 2000], and [Einstein, 2002] 
all require a numeric integration of splay, which is more costly to implement. 

The coefficient l/erf(^/2V5$-) is introduced into Eqs. (16 & 17) to force the 
trace of k (and therefore of K) to equal the trace of a 0 ® ao (viz., tr(ao <g> ao) = 1), 
which is the only non-zero invariant of tensor ao <8> ao- Alternatively, this coeffi- 
cient can be viewed as that scaling factor which is required to change the limits of 
integration from to in the Gaussian distributions present in these formulas. 


3.3.1 An Approximation The non-zero components of the analytic solutions to 
the intrinsic stiffness matrices listed in the appendices are expressed in terms of 
error functions with complex arguments (e.g., erf(z), z e C). Because this function 
is not found in most, standard, computer, math libraries, we introduce a simple ap- 
proximation to the analytic results derived in these appendices that one can readily 
employ; it being. 


K(g) * 


i(l + e" 2f2 ) 0 0 

0 |(l-e- 2f2 ) 0 

0 0 ^(l-e- 2 * 2 ) 


0</<l, 


(19) 


which is in keeping with the constraint that tr* = 1. Parameters / = 0 and 
/ = 1 apply to 2D splay with the normal to the membrane being in the 2- and 3- 
directions, respectively, while / = V 2 applies for 3D splay with transverse isotropy. 
Splay will be orthotropic whenever / £ V 2 ; specifically, there will be an elliptic 
symmetry in the transverse plane. 

The analytic (App. A) and approximate (Eq. 19) solutions for k are contrasted 
in Fig. 2 for 2D splay, while the formulae in App. B and Eq. (19) are contrasted in 
Fig. 3 for 3D splay with transverse isotropy. As a result, one can easily justify using 
the approximate solution stated in Eq. ( 19) over its analytic counterparts derived in 
the appendices, especially since selecting a Gaussian distribution to describe fiber 
dispersion was an assumption in the first place. 

Substituting Eq. (19) into Eq. (18) quantifies the global stiffness matrix K that 
appears in the elastic model of Eqs. (9 & 10) and its associated invariants in Eq. (6). 
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Figure 2 Plots of relative stiffness vs. standard deviation in the angle of fiber dispersion 
for 2D splay, as determined by Eq. (19) with / = 1 for the approximate solutions, and by 
Eqs. (A4 & A5) for the analytic solutions. 



Figure 3 Plots of relative stiffness vs. standard deviation in the angle of fiber dispersion for 
transversely isotropic 3D splay, as determined by Eq. (19) with / = y 2 for the approximate 
solutions, and by Eqs. (B4 & B5) for the analytic solutions. 
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4 Elasticity for Finite Elements 


Soft tissues are generally considered to be incompressible or nearly incompress- 
ible. In a finite element analysis that involves incompressibility, a standard dis- 
placement based interpolation method leads to ill-conditioning of the numerics 
[Malkus and Hughes, 1978]. Anticipating like difficulties, following [Flory, 1961], 
we seek to decouple pressure from displacement. This decoupling is compatible 
with two-field displacement-pressure interpolations that avoid volumetric locking 
in particular, and numerical ill-conditioning in general. 

Adopting the approach and notation of [Simo and Hughes, 1998, pp. 358-364], 
we define 


/ = det F = f , F = r' h F, C = F t F, B = FF T , (20) 


so that det F = 1, and therefore, det B = det C = 1. Consequently, the strain-energy 
density W is assumed to decouple as follows: 

6oW( C) = W(J) + W(C) + <W{ K, C), (21) 

where < W, <W, and *W are the dilational, distortional-isotropic, and distortional- 
anisotropic strain energies, respectively. 

The above definitions allow the general constitutive equation for elasticity 
stated in Eq. (8) to be recast as 


S = J 


30 


DEV 

'dW(C)' 

+ DEV 

d'Wi K, C)T 


3 C 


3 0 JJ 


so that, when pulled forward into the Eulerian frame, it becomes 
3W(0) , 


r = J 
wherein 


3 0 


dev 

[f d<W £) J 

+ dev 

F 0'iv(K,C) rx n 

l 

dC 


3 C JJ 


( 22 ) 


(23) 


DEV[*] = (•) - i tr((*)C)C 1 and dev[*] = (•) - | tr(*)l (24) 

are the respective Lagrangian and Eulerian deviatoric operators. Tensors F and 6>l 
are the deviatoric (volume preserving) and dilational (volume changing) parts of 
the deformation gradient F, respectively. Although 0 = J in a mathematical sense, 
we follow the admonition of Simo and Hughes and maintain their distinction, to 
remind ourselves that displacement and pressure are to be interpolated separately 
in finite-element codes suitable for soft-tissue analysis. 


4.1 A Simple Model 

The spherical strain-energy model advocated by [Simo and Hughes, 1998, pg. 361], 
and adopted, for example, by [Kaliske, 2000] , is 

<W{0) = k\{\{0 2 - 1) - In 0), (25) 

where k is the bulk modulus. Equation (25) is recommended because it is a convex 
function 4 , and because its gradient leads to a second-order accurate approxima- 

4 d 2 AVjdO 2 - |/c(l + 0 ~ 2 ) > 0, as © > 0 due to the conservation of mass 
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tion of the [Hencky, 1928] definition for dilatation 5 [Freed, 2004]. Equation (25) 
is appropriate for mean-dilation (or reduced/selective) integration schemes used in 
some finite-element codes. 

Equation (25) does not, however, fit the constraint criteria of the ujp finite- 
element scheme advocated by [Sussman and Bathe, 1987]. Thus, if this interpola- 
tion scheme is adopted, then 


'WO) = k\(J - l) 2 , (26) 

must be employed as the dilational strain-energy model. In these codes there is no 
distinction made between J and 0. 

Seeking an isotropic contribution to the deviatoric strain energy with attributes 
akin to those affiliated with the dilational part (viz., a convex function whose gra- 
dient produces a second-order accurate approximation of true strain), we assign 

2 < W / (C) = fi\(trC + trC -1 - 6), (27) 

where p is the shear modulus, with C -1 = {C 2 - tr(C)C + |[(tr C) 2 - tr C 2 ]l}/ detC 
from the Cayley-Hamilton theorem. Tensor C -1 exists because detC = 1; conse- 
quently, h{C) = ^[(trC) 2 - trC 2 ] = trC -1 = //(C -1 ). A [Mooney, 1940] material 
has different material constants assigned to invariants trC and trC -1 , in general, 
and in this sense, our model is a Mooney material of special form. 

For an anisotropic contribution to the deviatoric strain energy, going back to 
the precept that energy is the area under a force/displacement curve, we advocate 
that 

[tr(KC)]'/2 

^(K, C) = J o-(A) d A, (28) 

l 

where the fiber stress cr is allowed to be an arbitrary function of fiber stretch A; 
it is generally nonlinear in biological tissues. The upper limit of integration is the 
forth invariant, as it pertains to the deviatoric part of the deformed state. In order 
for this strain-energy function to be convex, it is necessary that E t (A ) > cr t (A ) for 
all A > 0, wherein 


dcr U) . o-(A) 

E t (A) = — — and cr t (A) - 


d A 


A 


(29) 


which are the fiber tangent-modulus and true-stress, respectively, with the fiber 
stretch A being quantified by 


A = (tr(KC))‘ /2 = (tr(FKF*))‘ /2 . 




(30) 


A physiologically based material model for cr(A) has recently been derived by 
[Freed and Doehring, 2004] that applies to crimped collagen fibers, which we have 
extended to meet our needs in App. C. 


5 e = lndet F ~ j(det F - det F J ) 
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Substituting the strain energies of Eqs. (25-28) into Eq. (22) produces an elas- 
tic constitutive model suitable for soft-tissue mechanics that when expressed in the 
Lagrangian frame becomes 

S = kJ \{0 - 6> -1 ) C" 1 + fiJ~ 2h DEV[±(I - C" 2 )] + o-,(A)J~ 2h DEV[K], (31) 

or equivalently, when substituted into Eq. (23), becomes 

T = kJ \{0 - 6T 1 ) I + n dev[i(B - B" 1 )] + cr t (A) dev[FKF T ], (32) 

when expressed in the Eulerian frame. The strain measure \{0 - 0~ l ) is a second- 
order accurate approximation of Hencky’s dilational strain field lndet F, while the 
strain tensor |(B - B" 1 ) is a second-order accurate approximation of the true-strain 
field In V [Freed, 2004]. Each contribution is convex provided that k > 0, jj. > 0, 
and E t (A) > cr t {A) for all A. 


4.2 Tangent Moduli 


The relationship between S and C in Eq. (31) is nonlinear. To obtain a finite- 
element solution with an iterative Newton-type solution process, that relationship 
must be linearized with respect to an incremental displacement. This involves the 
specification of a tangent modulus M. 6 To obtain this tangent, stress S is linearized 
over some time interval [t n , t n+ \] such that S„+i = S„ + M„dE with AE = E„ + i - 
E n = |(C„+i - C„) = F^E(x n ;t n ,r n+! )F n . Said differently, the tangent modulus 
corresponds to the slope affiliated with a forward-Euler integration step, and since 
it depends on step number n, it needs to be re-evaluated at each step along the 
solution path. 

Tensor M„ = 2dS/dC n = 4god 2 W/dC n dC n defines the tangent modulus in the 
Lagrangian frame, which can be pulled forward, component by component, into 
an Eulerian frame (commonly called the updated Lagrangian frame) according to 
the mapping m n ijkt = F n jj F n kK F n lL M n jJKL at the n* time step. Constructing the 
components of M„ in Voigt notation is addressed in App. D. 

From Eq. (21), suppressing the subscript n designating step number, the tan- 
gent modulus takes on the form 

M = M(<9) + M(C) + M(K, C), (33) 


where, through an application of the chain rule, one obtains 

- _ d 2 W(0) _ fd 2 'W(0) 80 80 8W(0) 8 2 0 \ 

~ 4 8C8C ~ { 80 2 8C ® 8C + 80 8C8C ) ’ 

- 8 2 <W( C) a ((8C^ 8 2 <W( C) 8C 5Ti>(C) d 2 C ' 

8C8C \\ 8C ) 8C8C 8C 8C 8080)’ 

5 2 TV(K, C) |7dC\ T . 8 2 <W( K,C) 8C f 8<W( K,C) ^ 2 C ) 

8C8C ~ [\^C/ ' 8080 ' 8C + 80 ' 8080)' 

(34) 


6 The tangent modulus M is often denoted as C, which we use to denote the right- 
deformation tensor of Cauchy. 
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with 


dCi 




dC kt 


and 


d 2 Cjj 

dCu 8C„ 


(35) 


= + i Q'O - i lk ij,c^ - (36) 


For the preferred volumetric strain energy in Eq. (25), which can be imple- 
mented into finite-element codes that use a reduced/selective integration scheme 
[Malkus and Hughes, 1978], one arrives at the following pressure tangent modu- 
lus 7 


M = 2 k 


. ? .<90 80 . i . 

0 + 0 } 


8 2 0 \ 
dCdCj' 


(37) 


where the gradients 80/8 C ® 80 /8C and 8 2 0/8C 8C are handled by the element 
technology of the particular finite element being employed. Alternatively, in an 
u/p formulation [Sussman and Bathe, 1987], the pressure tangent modulus corre- 
sponding to Eq. (26) is given by 8 


M = Kj(j C~ l ® + (1 - /)(2C _1 b C _1 - C _1 ® C -1 )), (38) 


where now the gradients 80/80 ® 80/80 and 8 2 0/8OdC are handled by us, the 
constitutive developers, as there is no distinction made between 0 and J in these 
codes. 

The isotropic tangent modulus corresponding to Eq. (27) is determined to be 


M = p{j 2h [0~ l ®0- 2 - ^C" 1 ® C~ 2 + C~ 2 ® C" 1 ) + ^(trC _1 )C _1 ® C _1 J 
+ \[{r lh trC - J 2h trC -1 )(C -1 H C _1 + ^C" 1 ® C" 1 ) 

- c 1 ® (,r 2/3 i - y 2/3 C" 2 ) - (r 2/3 i - f h c~ 2 ) ® cr 1 ]}. (39) 

Lastly, the anisotropic tangent modulus corresponding to Eq. (28) is found to 
be given by 


M = 7' 2/3 {(E f (T) - o- t (A))(r 2h r 2 K ® K 

- |(K ® C _1 + C _1 ® K) + \f h A 2 C" 1 ® C' 1 ) 

+ fcr,U)(/ /3 T 2 ( C^kiC" 1 + ^C _1 ® C -1 ) - K®C _1 -C _1 ® K)}. (40) 


7 These tensor gradients are useful when deriving the subsequent tangent moduli: 

dZ 8Z~ X , T 

— = I a I & — — = -Z aZ given that Z = Z 1 . 

oZ oZ 


8 The outer-dyadic tensor product = AijB ke , and the inner-dyadic tensor product 

|A h = \(A ik B j( + A i( B jk + A jk B i( + A jt B ik ), are established in App. D. 



16 


Alan D. Freed et al. 


5 Collagen Models 

The fiber tangent modulus E,(A) and true stress cr t (A) present in Eq. (40) are de- 
fined in Eq. (29), whose stretch A is established in Eq. (30). It is worth noting that 
the terms in the fourth-rank tensor representing the anisotropic tangent modulus 
are completely independent of the specific collagen fiber model chosen. This im- 
plies that a modification of fiber stress-strain law requires no re-formulation of the 
constitutive tangent matrix, making the model very flexible for different types of 
soft tissues. This is not the case with constitutive models represented by Eq. (1). 

In the Examples section that follows, two collagen stress-strain models are 
used. The first is the rule adopted by [Billiar and Sacks, 2000b] for aortic valve 
tissue; it being, 

cr(A) = A(e B{A2 ~ l)/2 - 1), (41) 

thus the two required terms in Eq. (40) are 

cr t (A) = AA _1 (e fi(A2 - I)/2 - 1) and E t {A) = A£Te S(/lM)/2 . (42) 

As an alternative to this phenomenological model, we also employ the structural 
model of [Freed and Doehring, 2004] that is based on the physiology of crimped 
collagen fibers. An adaptation of their algorithm is given in App. C (see Alg. 1). 
Specifically, given a fiber stretch A, this model returns the true stress cr /A and tan- 
gent modulus dcr/dT of the fiber. There are four physiologic parameters (material 
constants defined at the top of the algorithm) that the user must supply. 

To implement either of these two models in finite elements, or any other model 
for that matter, it is sufficient to change the lines of code corresponding to the 
scalar values E t (A) and cr t (A). This opens up the possibility of coding the material 
model once and selecting an appropriate fiber model with a passed parameter. 


6 Examples 

We offer our invariant theory as a computationally efficient alternative to statistical 
structural models of the type represented by Eq. (1). 

In the examples that follow, we match Eq. (31) to biaxial data for fresh aor- 
tic valve tissue. These data were generously supplied by Dr. Sacks from the En- 
gineered Tissue Mechanics Laboratory at the University of Pittsburgh, and have 
been reported on elsewhere [Billiar and Sacks, 2000a]. In the first example, we 
adopt an exponential fiber stress-strain rule — the same phenomenological model 
used by [Billiar and Sacks, 2000b]. In the second example, we adopt a structurally 
based collagen fiber model recently derived by [Freed and Doehring, 2004]. 

The original data were provided in the format of Lagrangian membrane ten- 
sion; that is, force per unit reference length. To characterize the constitutive model, 
it was necessary to convert the data to Lagrangian stress. Thus, a thickness of 
0.6 mm was assumed. It must also be noted that these data do not comprise a com- 
plete biaxial set, in the sense that there were finite off-axis deformation terms (i.e., 
F 12 t F 2 \ * 0, cf. Fig. 4) with no concomitant measurement of an off-axis stress. 
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Normalized x-coordinate (mm/mm) 

Figure 4 Mapping of a unit square to the current configuration, corresponding to the last 
data point of each of the seven protocols of [Billiar and Sacks, 2000a], Note that there exist 
finite F n + F 2 \ terms. 


As such, is was necessary to convert Eq. (31) into its first Piola-Kirchhoff stress 
(i.e. P = FS) counterpart for parameter estimation. 

An adaptive grid refinement (AGR) global optimization algorithm was imple- 
mented in Mathematica™ (Wolfram Research Incorporated, Champagne, IL), us- 
ing a commercially available global optimization algorithm (Global Optimization, 
Loehle Enterprises, Naperville, IL) [Doehring et al., 2004]. Parameters were si- 
multaneously fit to five separate biaxial load protocols, corresponding to fiber- 
to-cross-fiber membrane stress ratios of 30:60, 45:60, 60:60, 60:45, and 60:30 
N/m. These were the same protocols used to fit the data in the original publica- 
tions by [Billiar and Sacks, 2000a, Billiar and Sacks, 2000b]. During the estima- 
tion process, the membranes were considered to be incompressible and the shear 
modulus pertaining to the isotropic response was set to zero. The objective func- 
tion for the global minimum was defined as 



where P\\ corresponds to the fiber direction, P 2 i corresponds to the cross-fiber 
direction, and Nj represents the number of data points for the y* protocol. 
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6. 1 Exponential Fiber Model 


In this example, we specifically compare our constitutive model to both the data 
and to a statistical structural model with an exponential fiber stress that was pro- 
posed in [Billiar and Sacks, 2000b] and implemented by [Einstein, 2002] as 


S = i' 2/3 DEV 


J S f (d)R(6)N(d)®N(6)dO 


\~ n h 


(44) 


where the fiber stress-strain law S f{6) and the probability density function R(6) are 
given by 


S f (9) = - l) and R(0) = eXP (^W) (45) 

’ or V 2tt 

The fiber tangent modulus E t (A) and true stress cr t (A) for Eq. (31) are given by (42). 

Optimized parameters for the constitutive model based on our new invariant 
theory (Eqs. 31 & 42) were: A = 0.007, B = 21.6, and g — 0.795. Parameters for 
the constitutive model based on the statistical structure of Eqs. (44 & 45) were: 
A = 0.045, B = 21.0, and or = 0.192. Note that although the units for a in Eq. 
(45) are radians, our g parameter in Eq. (19) is not strictly an angle. To force the 
consistency constraint lim^ 0 Q k(s9Q t = ao<8iao required us to introduce a scaling 
factor of 1 /erf(/r/2 V2$-) into R(0). 

Overall, the fit was quite good for both models (see Fig. 5). The error calculated 
via Eq. (43) for the invariant model of Eqs. (31 & 42) was 0.40, as opposed to an 
error of 0.45 for the model of Eqs. (44 & 45). On average, Eqs. (31 & 42) tended 
to fit the fiber-direction stress slightly better, while Eqs. (44 & 45) tended to fit 
the cross-fiber stress slightly better. With regard to our constitutive model only, the 
toe-region of the stress- strain curve (commonly viewed as a transition between the 
extinction of collagen crimp and the linear behavior of straightened collagen) was 
slightly under-predicted in the fiber direction and over-predicted in the cross-fiber 
direction. 


6.2 Crimped Collagen Model 

Because the modulus E t (A) and true stress ay (A) are generic scalar components 
of the anisotropic tangent modulus in Eq. (40), we are free to adopt any reason- 
able fiber stress-strain rule, without making the complexity due to tangent modulus 
construction prohibitive. This is a desirable precondition for finite-element analy- 
sis. 

As our second example, consider the micro- structurally based collagen fiber 
stress-strain rule based on the physiology of crimped collagen fibers that was re- 
cently proposed by [Freed and Doehring, 2004] (Alg. 1 in App. C). This is an algo- 
rithm for the elastic response of crimped collagen fibers, based on the observation 
that fibril crimp has a three-dimensional structure at the pm scale whose geometry 
can be approximated as a cylindrical helix. For pre-failure analysis, the model is 
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Deformation (FI 1,F22) Deformation (F 1 1 ,F22) 




1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 

Deformation (F 1 1 ,F22) 


Figure 5 Constitutive model fit to biaxial data from fresh aortic valve tissue: exponential 
fiber stress-strain rule (left: fiber direction, right: cross-fiber direction). Dots are the original 
data. Solid lines are from our constitutive model Eqs. (31 & 42). Dashed lines are from the 
statistical model in Eq. (44) with Eq. (45) for the stress/stretch law. 
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defined in terms of three physiologic parameters: the initial normalized wavelength 
of the crimp H 0 /r 0y the initial normalized amplitude of the crimp R 0 /r 0 , and the 
elastic modulus of the collagen fiber in the linear region E f. Parameters Hq and Rq 
are normalized with respect to fibril radius r 0 . 

Values for these parameters were estimated at: H 0 /r 0 = 14.4, Ro/r 0 = 2.19, 
Ef = 10.6MPa, and g = 0.768. Overall, the fit of this variation of our constitutive 
model was excellent (see Fig. 6), with an error of 0.37. Only the predictions for 
the sixth protocol have a significant difference from the data. 


7 Finite Element Implementation 

Equation (31) and Alg. 1 from App. C was implemented into Adina™ (Adina 
R&D Inc., Watertown, MA). A lOx 10 mm square of tissue was meshed with 1764 
nodes and 1200 linear solid elements with constant pressure interpolation. A two- 
field pressure/displacement interpolation was utilized [Sussman and Bathe, 1987]. 
Displacement boundary conditions were applied to edge nodes. To characterize the 
full range of behavior, the load was divided into 100 load steps. Simulations were 
performed for three biaxial stretch ratios: X \\ : A .22 = 1.2 : 1.4, 1.1 : 1.4, and 1.0: 1.4. 
The solutions were obtained using a full-Newton method with a sparse equation 
solver. Cpu times were 154, 149, and 150 seconds run on a Linux box containing 
a single 2.4 GHz Pentium IV processor. In addition to the parametric values stated 
earlier, the following isotropic moduli were assigned: k = 20MPa and [i = 500Pa. 
There was excellent agreement between the FEA and theoretical solutions (see 
Fig. 7), with essentially no detectable error. 


8 Conclusion 

We have proposed an efficient, invariant-based alternative to structural constitutive 
equations that accounts for a statistical dispersion of fibers. In contrast to existing 
models, our new invariant theory easily handles a 3D fiber population with a sin- 
gle mean preferred direction. The invariant theory is based on a novel closed-form 
‘splay invariant’ that requires a single parameter in the 2D case, and two parame- 
ters in the 3D case. The model is polyconvex, and fits biaxial data for aortic- valve 
tissue better than existing aortic- valve models. A modification in the fiber stress- 
strain law requires no re-formulation of the constitutive tangent matrix, making 
the model flexible for different types of soft tissues. Most importantly, the model 
is computationally expedient in a finite-element analysis. 

Acknowledgements The authors take this opportunity to thank Prof. Michael Sacks at 
the University of Pittsburgh for providing us with his experimental data, and to Dr. Todd 
Doehring and Mr. Dimitri Deserranno at the Cleveland Clinic for many delightful discus- 
sions on this and related topics. 
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Deformation (FI 1 ,F22) Deformation (FI 1 ,F22) 



I 1.1 1-2 1.3 1.4 1.5 1.6 1.7 1.8 1 1.1 1.2 1.3 1.4 1,5 1.6 1.7 1,8 

Deformation (F11.F22) Deformation (FI 1,F22) 



1 1.1 1-2 1.3 1.4 1.5 1.6 1.7 1.8 

Deformation (FI I, F22) 


Figure 6 Constitutive model fit to biaxial data from fresh aortic-valve tissue: fiber crimp 
stress-strain rule (left: fiber direction, right: cross-fiber direction). Dots are the original data. 
Solid lines are our constitutive model Eq. (31) using Alg. 1 from App. C for the stress/stretch 
law. 
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Figure 7 Finite-element simulation of three biaxial stretch ratios: A n : A 22 = 

1.2: 1.4, 1.1 .1.4, and 1.0: 1.4 (solid lines). Left: fiber direction, right: cross-fiber direction). 
Dots are theoretical solutions. 

Appendices 

The primary reason for adopting a Gaussian distribution to describe fiber splay is 
that the corresponding stiffness matrix k can be computed analytically. Alterna- 
tively, [Hurschler et al., 1997] have employed a von Mises distribution for splay 
in conjunction with a Weibull distribution for crimp that they collectively solve 
numerically. 


A Two-Dimensional Splay 


We recall that our local coordinate system was chosen so that 


a 0 



( cos 

tf = < sin# l , 

1 0 J 


(Al) 


and as such 


e /® e / 


cos 2 0 sin 8 cos 8 0 
sin 6 cos 8 sin 2 8 0 , 

0 0 0 


(A2) 


which leads to an expression for Eq. (16) that can be solved analytically; it being, 


Kn 0 0 
0 K22 0 , 
0 0 0 


K{g) = 


(A3) 
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where the two non- zero stiffness components have values of 


*" = 2 + 


2 + 


= 2 " 




wherein i = V-I is the unit imaginary number. 

The sum of two error functions whose arguments are complex conjugates, i.e., 
erfU + iy ) + erf(* - iy), produces a real result, and as such, Ku and K 22 are both 
reals. 


B Three-Dimensional Splay 

Here the formulation is somewhat different; specifically, 

fl'j ( cos 9 \ 

ao = io>, e/ = jsin# cos^i , (Bl) 

[oj (sin# sin 0 ] 

and as such 

cos 2 # sin # cos 6 cos 0 sin 8 cos 6 sin 0 
e/®e/= sin 9 cos 6 cos 0 sin 2 9 cos 2 0 sin 2 9 sin 0 cos 0 , (B2) 

sin# cos# sin 0 sin 2 # sin 0 cos 0 sin 2 # sin 2 0 

which leads to an expression for Eq. (17) that can be solved analytically; it being, 

k u 0 0 

*(?) = 0 K 22 0 , (B3) 

. 0 0 /C33 

where the two non-zero stiffness components have values of 


*>> = 2 + 


M*) 


which is the same as Eq. (A4), and 


*22 = *33 = 7 

4 


" 4 8erf te) 


where this value for K 22 is exactly half that of Eq. (A5), wherein K 33 = 0, but here 

*33 = * 22 - 
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Algorithm 1 

Given Ho/ro, Rq/ r 0 , Ef, and A u , where Hq is the initial wavelength of crimp, 

Rq is the initial amplitude of crimp, r 0 is the initial fibril radius, E f is the elastic 
modulus of the fiber in the linear region, and A u is its ultimate stretch, then: 

Set r 0 = 1 so that H 0 = H 0 /r 0 & R 0 = R 0 /r 0 . 

Compute the constant parameters: 

Lo = ((2nR 0 ) 2 + H 2 ) 1 ' 2 , A = Lq/Hq < A u , 

E s = E f /{(H 0 /Lo) 2 + [H 0 (Lo - H 0 )/L 2 ][ 1 + 37/6^ + 2(Lo/nr 0 ) 2 ]}, 
where Lq is the chord length of helix over one wavelength, while A is the stretch 
and E s is the secant modulus at the transition between the toe and linear regions. 
If A < A Then 

£ = 6(7tr 0 ) 2 \A 2 + (4/ r 2 - l)d 2 ]d 

/ {A{3H 2 (A 2 - A 2 )[3A 2 + (87T 2 - 3 )A 2 ] + 8(^r 0 ) 2 [10A 2 + (3 tt 2 - 10)d 2 ]}}, 
df/d A = {\SH 2 0 (nr 0 ) 2 [3A 6 + (28 n 2 - 3)A 4 A 2 + (32tt 4 - 8^ - 3)A 2 d 4 
+ (32^ 4 - 20n 2 + 3)4 6 ] + 48(^r 0 ) 4 [10A 4 
+ (1 17 tt 2 - 20) A 2 d 2 + (12^ 4 - 43 ^ 2 + 10)d 4 ]} 

/ [A{3H 2 (A 2 - A 2 )[3A 2 + (8 n 2 - 3)d 2 ] 

+ 8(^t 0 ) 2 [10A 2 + (3; r 2 - 10)d 2 ]} 2 }, 
a/ A = E£{A - 1)/A 2 & dcr/dA = E s [£/A 2 + (df/dA)(A-l)/A] 

Else If A < A < A u Then 

<rlA = E s (A-\)IAA + E f {A-A)IA & dcr/dd = E f 
Else Fibril Failure 

cr/A = 0 & d cr/dA = 0. 

Return cr/ A and dcr/dd. 


C Collagen Model 

A micro-structural model based on the physiology of crimped collagen fibers was 
recently derived by [Freed and Doehring, 2004] that we have altered to meet our 
needs (see Alg. 1); specifically, given a fiber stretch A, this model returns the true 
stress c/A and tangent modulus dcr/ dA of the fiber. There are four physiologic 
parameters (material constants defined at the top of the algorithm) that the user 
must supply; three if failure is not to be considered. 

A typical pair of response curves are plotted in Fig. 8. The discontinuity ob- 
served in the do-jdA curve indicates that the model presented in Alg. 1 predicts a 
stress cr response that is continuous and once-differentiable in stretch A up to fiber 
failure, which is not depicted in this figure. The values assigned to the model to 
obtain these curves were Ho/ r 0 = 27.5, R 0 /r 0 = 2, and Ef = 45MPa, which results 
in a transitional stretch of A « 1.1. 

This is but one example of a fiber stress/stretch model. 


D Voigt Notation 


Because the Lagrangian fields for stress and strain £)y, i, j = 1, 2, 3, are sym- 
metric tensors, it is customary to express their components as six-dimensional ar- 
rays S a and E (X , a = 1, 2, . . . , 6. These arrays are not vector fields in the sense that 
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X 

Figure 8 Typical soft-tissue response curves for cr/A and dcr/dA. 


they do not obey the tensor-transformation law. Nevertheless, they have proven to 
be very useful in the employ of finite elements. 

Given a Cartesian reference frame, the stress and strain arrays can be expressed 
in terms of their tensor components via 9 



Sn 

S22 

S 33 

512 = S21 

513 = S 31 

S23 = S 32 


and 


'EA ( E 11 ' 

E 2 E 22 

£3 £33 

P — ) ' K ~ \ 

“ £4 E n + E 2 i ’ 

£5 £13 + £31 

E(, E22 + £32, 


(Dl) 


which is commonly referred to as Voigt notation [Belytschko et al., 2000, pp. 615- 
618]. From thermodynamics, stress is described in terms of a potential function W 
through the gradient 


Sij = go 


I dW cW\ 

Uc, 7 + dcj 


Sij = Sji, (D2) 


where C,- ; is the Lagrangian (or Green) deformation field. It is more common to see 
the above formula written in the condensed notation of tensors S = 2godWjdO (see 
Eq. 8 ) with the implication being that its components are computed via Eq. (D2). 

9 A note of caution. Some implementations exchange S 4 and S 6 , and likewise, £4 and E 6 . 
This has no effect on formulae written using Voigt notation, provided that the components 
are correctly mapped between their tensor and Voigt representations. 
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A stress increment d S,j is therefore related to a strain increment d E tj = ±d Qj 
in the Lagrangian frame through the linear approximation 


d Sij = M ijke dE ke with 


Miju - 


dSjj 

dE u 


d S a = M a pdEp, (D3) 


where the tangent modulus Mijm has a Voigt representation of 


M\m Ml 122 Mii33 Mm2 Mm3 Ml 123 

M2211 M2222 M2233 M2212 M2213 M2223 
M3311 M3322 M3333 M3312 M3313 M3323 
M1211 M1222 M]233 M1212 M1213 M1223 
M1311 Mi 322 M1333 M1312 Mj3i3 Mi 323 
M2311 M2322 M2333 M2312 M2313 M2323 

(D4) 

with M a p t Mp a unless the Miju possess major symmetry M ijki - M k uj. The tan- 
gent modulus always possesses minor symmetries M, ; h = M i; ^ = Mjm = Mju k 
because of the inherent symmetries in S t j and C,y. Major symmetry follows auto- 
matically if stress Sij is given by a potential function F such that Sij = dFjdEij, as 
is the case in elasticity (see Eq. D2). 

From Eqs. (D2 & D3), it follows that the components of M^e are obtained via 


M\\ M \ 2 iK/13 M\4 M 15 M 16 


M21 M22 M23 M24 M25 M26 


M31 M32 M33 M34 M35 M36 


M41 M42 M43 M44 M45 M46 


M51 M52 M53 M54 M55 M56 


_Mgi Mfi 2 M<53 Mg 4 M65 Mgg 



/ d 2 W d 2 W d 2 W d 2 W \ 

[ dCijdCu + dC u dCa + dCjidCu + dCjidC ek ) ’ 


(D5) 


so that M l]kt = Mijtk = Mji k ( ~ M j ikk and M, ; k = M k aj. In tensor notation, the 
tangent modulus is usually defined as M = 4g 0 d 2 W/dCdC (see Eq. 34) with the 
implication being that its components are computed via Eq. (D5). 


D.l Tensor Products 


Two tensor products naturally arise in constitutive construction that we denote as 
A B B and A® B, and call the inner- and outer-dyadic products, respectively, or the 
circle- and box-products for short, wherein A and B are taken to be symmetric. 
The box product A E B is defined by the sum 


(A B B)ij k c - j ( Ai k B j( + Ai(B j k -1- A j k B ik + A j{B ik ), (D6) 


which possesses both minor and major symmetries, therefore (A&B) a p = (A^B)p a , 
and consequently, this product is commutative, i.e., AaB = B E A. The Voigt 
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representation of (A E B) a p has the symmetric matrix components: 

r (A E fi)n = A\\B\\, 

(A El B) 12 = A\iB\2, 

(A E B) 13 = A13B13, 

(A E 5)14 = |(Anfii2 + A12B11), 

(A E 5)15 = ^(Aiifii 3 + A U B U ), 

(A E B) 16 = ^(Ai 2 fil 3 + A13B12), 

(A E B ) 22 = A22B22, 

(A E B ) 23 = A23B23, 

(A E 5)24 = 5 (A 12^22 + A22B12), 

(A E B )25 = ^ (A 1 2 B23 + A23512), 

< (A E 5)26 = 5 ^ 22 ^ 23 + A23B22), (D 7 ) 

(A E 5)33 = A33B33, 

(A E 5)34 = ^ (A 1 3 B23 + A23513), 

(A E 5)35 = ^ (A 13 533 + A33513), 

(A E B) 36 = ^ (A 23 533 + A33B23), 

(A E 5)44 = ^(An5 2 2 + 2A12B12 + ^22^11), 

(A E 5)45 = ^ (A 1 1 Z ?23 + ^ 13^12 + A\2B\3 + A23B\\), 

(A E 5)46 = ^(Ai 2523 + ^ 13^22 + ^ 22^13 + ^ 23 ^ 12 ), 

(A E 5)55 = \{A\\B 33 + 2 A 13 5 i 3 + A23B11), 

(A E B) 56 = ^ (A 12^33 + Ai 3 i ?23 + ^ 23^13 + ^ 33 ^ 12 ), 

„ (A ® 5) 6 6 = 5(^22533 + 2 A 23^23 + A 33 522 )' 

The circle product A 0 B has components 

(A 0 B)iju = j(A (J - 5 jtr + A^B^ + A + A ;i - 5 ^) = A^B^e, (D8) 

which reduce down to the simple expression (A 0 6), 7 k = AijBu as a consequence 
of A and B being symmetric. This well-known product has a Voigt notation of 


Aii5n A n 5 2 2 Ah 5 3 3 A n 5i2 A\\B\3 Aii#23 


(A 0 B) ap = 


A 22^11 A22B22 A22B33 A22B12 A 22^13 A 22^23 
A335ii A33B22 A33B33 A32B12 A32B13 A33B23 

AnBn A12B22 A12B23 A n B n Ai2#i3 A\2B23 

Ai35n Ai3#22 A13533 A13B12 Ai 3 5i3 Ai3#23 
A235 h A23B22 A23B33 A23B12 A23B13 A23B23. 


(D9) 


where (A 0 B ) Qjg is not symmetric unless A = B; however, (A 0 B) + (B 0 A) 
does yield a symmetric matrix in Voigt notation. This is because A 0 B possesses 
minor symmetry, but not major symmetry, in general. This dyadic product is not 
commutative. 
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D.2 Example 

The theory of linear elasticity has stress components 

S\ = Sn = T(£n + £22 + £33) + 2 /i£n, 

5 2 = S22 = /i(£ii + £22 + £33) + 2 pE 22 , 

53 — 533 = T(£n + £22 + £33) + 2 a*£ 33 , 

5 4 = S 12 = 2pEi2, 

55 = 5n = 2>u£ 13, 

,56 = 523 = 2pE 2 3, 

and as such, its tangent modulus is readily determined to be 
M a p = A{1 ® I) ap + 2/^(7 B 7) a/J , 
wherein A and p are the elastic moduli, and where 
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1 

100 

0 
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0 

0 

000 
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0. 


0 

00 
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V 2 . 


The 1 / 2 ’s in (/ b I) a/3 are offset by the 2’s present in d Ep. 
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